subroutine GL(ss,x,y,energy_out)
use var
use type
use parallel
implicit none
include 'mpif.h'
integer :: i,ii,x,y,iout, ss, tt
real(dp) :: delta, delta2, dmc_cur, temp_cur
real(dp) :: energy(5), E1, E2, DE, ac, dt_cur, max_ac, min_E
real(dp), intent(out) :: energy_out(5)
real(dp),dimension(:), allocatable :: cur_energy
logical :: flag = .true.
if (ID == Master) then
   allocate(cur_energy(5*size))
endif
iout = 300 + ID + 1

delta = (dt(2)-dt(1))/dt(3)
if (x/=0) then
   call free_energy(energy)
   write(*,'(I5,I8,2X,4(f15.6),2(f20.6))')ID+1,x,energy,Escale*energy(5)/kbT
   dt_cur = dt(1) + delta*y
   if (y>dt(3)) dt_cur = dt(2)
   call  ginzburg_landau(dt_cur)
   x = x + 1
   y = y + 1
endif
E1 = 0   ! Initial Energy
ac = 1
max_ac = 1
flag = .true.
tt = 1
Do while (flag)
   if (mod(x,sample)==0) then
      call free_energy(energy)
      E2 = energy(5)
      DE = (E1-E2)
      ac = abs(DE/E1)*100
      E1 = E2
      write(iout,'(I8,2X,4(f15.6),2(f20.6))')x,energy,Escale*energy(5)/kbT
      CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
      CALL MPI_ALLREDUCE (ac, max_ac, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, IERR)
      CALL MPI_GATHER(energy,5,MPI_DOUBLE_PRECISION,cur_energy,5,MPI_DOUBLE_PRECISION,master,MPI_COMM_WORLD,ierr)
      if (ID == Master) then
         open (22, file = 'timestep')
         Do ii=1,size
            write(22,'(I5,I8,2X,4(f15.6),2(f20.6))')ii,x,cur_energy(1+(ii-1)*5:ii*5),Escale*cur_energy(ii*5)/kbT
            write(*,'(I5,I8,2X,4(f15.6),2(f20.6))')ii,x,cur_energy(1+(ii-1)*5:ii*5),Escale*cur_energy(ii*5)/kbT
         enddo
         close(22)
      endif
      trj_E(:,tt) = energy
      tt = tt+1
      if (max_ac < accuracy(2)) flag = .false.

   endif
   if (mod(x,traj(1))==0) then
      if (ID == Master) then
         Do ii=1,size
            write(*,'(I5,I8,2X,4(f15.6),2(f20.6))')ii,x,cur_energy(1+(ii-1)*5:ii*5),Escale*cur_energy(ii*5)/kbT
            write(222,'(I5,I8,2X,4(f15.6),2(f20.6))')ii,x,cur_energy(1+(ii-1)*5:ii*5),Escale*cur_energy(ii*5)/kbT
         enddo
      endif
      if (traj(2)==1) call wtraj(ss,x,y)
!       call qmga(0,x)
      tt = 1
   endif
   
   dt_cur = dt(1) + delta*y
   if (y>dt(3)) dt_cur = dt(2)
   call  ginzburg_landau(dt_cur)
   
   x = x + 1
   y = y + 1
enddo
if (ID == Master) then
Do ii=1,size
   write(222,'(I5,I8,2X,4(f15.6),2(f20.6))')ii,x,cur_energy(1+(ii-1)*5:ii*5),Escale*cur_energy(ii*5)/kbT
enddo
endif
energy_out = energy
endsubroutine GL
